#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <TCanvas.h>
#include <TF1.h>
#include <TFile.h>
#include <TGraph.h>
#include <TH2F.h>
#include <TKey.h>
#include <TMath.h>
#include <TTree.h>
#include <TVector3.h>

void ge_dc_vs_cos_theta_plot(TCanvas *, std::string const &);

void ge_dc_vs_cos_theta_plot(TCanvas *canvas, std::string const &filename)
{
  TFile file(filename.c_str());
  if (!file.IsOpen()) {
    throw std::runtime_error("Could not open file.");
  }
  TObject *obj = file.GetKey("is430")->ReadObj();
  if (!obj || 0 != strcmp("TTree", obj->ClassName())) {
    throw std::runtime_error("Missing tree 'is430'.");
  }
  TTree *tree = dynamic_cast<TTree *>(obj);

  float ge_e;
  UShort_t ge_id;
  float ge_x;
  float ge_y;
  float ge_z;
  tree->SetBranchAddress("ge_e", &ge_e);
  tree->SetBranchAddress("ge_id", &ge_id);
  tree->SetBranchAddress("ge_x", &ge_x);
  tree->SetBranchAddress("ge_y", &ge_y);
  tree->SetBranchAddress("ge_z", &ge_z);

  UChar_t fb_channel;
  float fb_gamma, fb_beta;
  float fb_x, fb_y, fb_z;
  tree->SetBranchAddress("fb_channel", &fb_channel);
  tree->SetBranchAddress("fb_beta", &fb_beta);
  tree->SetBranchAddress("fb_gamma", &fb_gamma);
  tree->SetBranchAddress("fb_x", &fb_x);
  tree->SetBranchAddress("fb_y", &fb_y);
  tree->SetBranchAddress("fb_z", &fb_z);

  TH2F hist("hist", "ge_vs_ne_theta", 20, -1, 1, 8e3, 0, 8e3);

  float gamma = 2.85 / 939.0 + 1;
  float beta = sqrt(1 - 1 / (gamma * gamma));
  for (int i = 0; tree->GetEntries() > i; ++i) {
    tree->GetEntry(i);
//    float cos_theta = gx * fb_x + gy * fb_y + gz * fb_z;
    float cos_theta = ge_z;
    float E0 = ge_e * gamma * (1 - beta * cos_theta);
cos_theta = ge_x * fb_x + ge_y * fb_y + ge_z * fb_z;
    hist.Fill(cos_theta, E0);
  }
  canvas->cd(0);
  hist.DrawCopy("zcol");
}
